Phase separation in the crust of accreting neutron stars 
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Nucleosynthesis, on the surface of accreting neutron stars, produces a range of chemical elements. 
We perform molecular dynamics simulations of crystallization to see how this complex composition 
forms new neutron star crust. We find chemical separation, with the liquid ocean phase greatly 
enriched in low atomic number elements compared to the solid crust. This phase separation should 
change many crust properties such as the thermal conductivity and shear modulus. 
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I. INTRODUCTION 

Phase separation is important for white dwarf stars 
As a star cools, and crystallization takes place, the 
crystal phase is enriched in oxygen while the liquid is 
enriched in carbon. We believe phase separation may 
also be important for neutron stars that accrete material 
from a companion. This material can undergo nuclear 
reactions involving rapid proton capture (the rp process) 
to synthesize a variety of medium mass nuclei |2f . Fur- 
ther accretion increases the density of a fluid element 
until crystallization occurs. However, as we explicitly 
demonstrate with molecular dynamics simulations, crys- 
tallization is accompanied by chemical separation. The 
composition of the new solid crust is very different from 
the remaining liquid ocean. This changes many proper- 
ties of the crust and can impact many observables. 

With chemical separation, the liquid ocean, see Fig. 
[TJ is greatly enriched in low atomic number Z elements. 
Carbon, if present, may be depleted in the crystal (crust) 
and enriched in the liquid ocean phase. Also, chemi- 
cal separation may change the thermal conductivity of 
the crust and its temperature profile. Indeed some neu- 
tron stars are observed to produce energetic X-ray bursts 
known as superbursts. These are thought to involve the 
unstable thermonuclear burning of carbon 0, 0, • How- 
ever it is unclear how the initial carbon concentration is 
obtained and how the ignition temperature is reached. 

Chemical separation may significantly change the 
thickness, shear modulus, and breaking strain of the 
crust. This could change the shape of a neutron star and 



the radiation of periodic gravitational waves @ pj . Fur- 
thermore, changes in the crust could change the proper- 
ties of quasi periodic oscillations that may be observable 
in thermonuclear bursts. 
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FIG. 1: (Color on line) Schematic diagram of the surface of 
an accreting neutron star. This paper focuses on chemical 
separation upon crystallization at the boundary between the 
liquid ocean and the solid crust. We find that the ocean is 
enriched in low Z elements. Note that the boundary between 
the inner and outer crust is not shown. 



The process we consider here is distinct from sedimen- 
tation, which may also occur at lower densities where 
the ions are fluid Q. However, we are not aware of any 
previous calculations of chemical separation from crystal- 
lization for accreting neutron stars. Jones has considered 
how a range of compositions may change the properties 
of the crust of a non-accreting neutron star. In ref. Q he 
considers phase separation, but only based on quite early 
work on the free energy of the two-component Coulomb 
plasma. Instead, Jones suggests that the system will form 
an amorphous solid [Iol |. 
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The ash resulting from rapid proton capture (rp pro- 
cess) nuclear reactions is expected to have a complex 
composition involving a number of different chemical el- 
ements @, [ll[ . Unfortunately, it can be difficult to con- 
struct the phase diagram for a multi-component system. 
The pure one component plasma (OCP) phase diagram 
is well known. The liquid solidifies when the ratio of 
a typical Coulomb energy to the thermal energy kT is 
r « 175 p|. The parameter T is defined, 
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where the ion charge is Ze, the temperature is T, and the 
ion sphere radius a describes a typical distance between 
ions, a = (3/47rn) 1/3 . Here n is the ion (number) density. 
The phase diagram for binary mixtures has also been 
determined. See for example [15| . For binary mixtures, 
the solid phase is enriched in the high Z ion and the 
liquid phase is enriched in the low Z ion. 

Often, the theoretical phase diagram is constructed 
from extremely accurate calculations of the free ener- 
gies of the solid and liquid phases. The melting point 
is determined by equating these two free energies. Very 
accurate calculations are needed because the free ener- 
gies are nearly parallel as a function of T. A small error 
in the free energy of one phase can lead to a large error 
in the melting point. Therefore, it may be very difficult 
to compute the free energy of multi-component systems 
with enough accuracy to determine the phase diagram. 

Instead, in this paper we determine the phase diagram 
of a multi-component system directly via molecular dy- 
namics (MD) simulations. Our simulation volume con- 
tains regions of both the liquid and solid phase. This 
approach has many advantages. It is simple and robust. 
Delicate free energy calculations are not needed. One can 
directly measure the composition of the two phases that 
are in equilibrium. Furthermore, one can run simulations 
with arbitrarily complicated compositions. 

However, there are two limitations to direct molecular 
dynamics simulation. First, finite size effects may be sig- 
nificant because a large fraction of the ions are near the 
interfaces between the two phases. We minimize finite 
size effects by using a moderately large number of ions 
27,648 and we measure the composition of the two phases 
in regions that are away from the interfaces. Second, it 
can take a long time for the two phase system to come 
into thermodynamic equilibrium. We address non equi- 
librium effects by running for a total simulation time of 
151 million fm/c (over six million MD time steps) and 
by monitoring the time dependence of the composition 
of the two phases. Still, as we discuss below, the sys- 
tem may not be in full equilibrium and this may be an 
important question for further work. Nevertheless, we 
start with equal compositions and find dramatically dif- 
ferent compositions for the liquid and solid, where the 
difference has increased with simulation time. If the sys- 
tem is not in equilibrium by the end of our simulation, 
we expect the difference between liquid and solid to only 



increase further with time. Therefore, we do not think 
non-equilibrium effects will change our conclusion that 
the liquid and solid have very different compositions. 

This paper is organized as follows. In Section [H] we 
describe our molecular dynamics simulation. Results are 
presented in Section IIIII and we conclude in Section IIV1 



II. MOLECULAR DYNAMICS SIMULATION 

We now describe the initial composition for our sim- 
ulation. Schatz et al. have calculated the rapid proton 
capture (rp) process of hydrogen burning on the surface 
of an accreting neutron star [2| . This produces a variety 
of nuclei up to atomic masses A « 100. Gupta et al. 
p| then calculate how the composition of the rp process 
ash evolves because of electron capture and light particle 
reactions as the material is buried by further accretion. 
Their final composition, at a density of 2.16 x 10 11 g/cm 3 
(near neutron drip at the bottom of the outer crust) has 
forty percent of the ions with atomic number Z = 34, 
while an additional 10 % have Z = 33. The remaining 
half of the ions have a range of lower Z from Z = 8 to 
32. Finally there is a small abundance of Z = 36 and 
Z = 47. 
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FIG. 2: (Color on line) Abundance (by number) of chemical 
elements versus atomic number Z. The plus symbols show the 
initial composition of the mixture. The final compositions of 
the liquid phase, open green circles, and solid phase, filled red 
squares, are show after a simulation time of 151 x 10 6 fm/c, 
see Section 11111 



For simplicity we use the Gupta et al. abundances be- 
cause we have them available. However these abundances 
were calculated assuming no phase separation. Therefore 
they have not been determined self-consistently if there is 
phase separation. Nevertheless, we use them to provide 
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TABLE I: Abundance y z (by number) of chemical element Z. 
Results are presented for the original mixture and for the final 
liquid and solid phases after a simulation time of 151 x 10 6 
fm/c, see text. 



z 


Mixture 


Liquid 


Solid 


8 





.0301 





.0529 





.0087 


10 





.0116 





.0205 





.0021 


12 





.0023 





.0043 





.0006 


14 





.0023 





.0043 





.0005 


15 





.0023 





.0043 





.0004 


20 





.0046 





.0055 





.0029 


22 





.0810 





.1024 





.0616 


24 





.0718 





.0816 





.0635 


26 





.1019 





.1065 





.1017 


27 





.0023 





.0025 





.0027 


28 





.0764 





.0744 





.0746 


30 





.0856 





.0773 





.0949 


32 





.0116 





.0099 





.0130 


33 





.1250 





.1079 





.1388 


34 





.3866 





.3408 





.4297 


36 





.0023 





.0021 





.0030 


47 





.0023 





.0030 





.0013 



a first orientation. Note that we use abundances calcu- 
lated near 10 11 g/cm 3 , while the ocean/ crust boundary 
may be near 10 10 g/cm 3 . The differences in composition 
at these two densities may be primarily do to a modest 
amount of electron capture. This should not significantly 
change our results. Perhaps phase separation will lead 
to more important changes in the abundances. Chemi- 
cal separation is expected to change compositions over a 
large range of densities in addition to densities near the 
ocean/ crust interface. For example, changes in compo- 
sition of the liquid, near the crust interface, are expected 
to diffuse throughout the ocean. As we discuss in section 
HVI futurc calculations of abundances including phase sep- 
aration would be very useful. 

As an initial composition we chose 432 ions with Z and 
mass number A drawn at random according to the Gupta 
et al. abundances. This is shown in Fig. fj] and listed in 
Tableland closely approximates the original distribution 
up to limitations of the small statistics. We chose such a 
small system, 432 ions, to simplify producing the original 
solid configuration, see below. Note that the liquid and 
solid phase results shown in Fig. [2] will be discussed in 
Section [ml 

At these densities, electrons form a relativistic degen- 
erate Fermi gas. The ions are fully pressure ionized and 
interact with each other via screened Coulomb interac- 
tions. The potential between the ith and jth ion is as- 
sumed to be, 

Vij{r) = ^M e -r/A, (2) 

Here the ion charges are Zi and Zj, r is their separation 
and the electron screening length is A. For cold rela- 
tivistic electrons, the Thomas Fermi screening length is 
A -1 = 2a 1 l 2 kp i 'it 1 / 2 where the electron Fermi momen- 
tum kp is kp — (37r 2 n e ) 1 / 3 and a is the fine structure 



constant. Finally the electron density n e is equal to the 
ion charge density, n e — (Z)n, where n is the ion density 
and (Z) is the average charge. Note that we are inter- 
ested in temperatures near the melting point where the 
ion thermal de Broglie wave length is much shorter than 
the inter ion spacing. Therefore quantum corrections to 
the ion motion should be very small. 

We now describe the initial conditions for our classical 
MD simulation. It can be difficult to obtain an equilib- 
rium crystal configuration for a large system involving a 
mixture of ions. Therefore, we start with a very small 
system of 432 ions with random coordinates at a high 
temperature and cool the system a number of times by 
re-scaling the velocities until the system solidifies. Here 
the velocities of all of the ions are multiplied by a com- 
mon factor so that the kinetic energy per ion is 3T/2 for 
a series of decreasing temperatures T. Next four copies 
of this solid configuration were placed in the top half of a 
larger simulation volume along with four copies of a 432 
ion liquid configuration. The resulting system with 3456 
ions was evolved in time until it fully crystalizcd. Finally, 
four copies of this 3456 ion crystal were placed in the top 
half of the final simulation volume along with four copies 
of a 3456 ion liquid configuration. This final system has 
27648 ions and consists of a solid phase above a liquid 
phase. Note that the initial compositions of these two 
phases are equal. 

Our results can be scaled to different densities. For 
historical reasons, our simulation was run at a relatively 
high ion density of n = 7.18 x 10~ 5 fm -3 . This corre- 
sponds to a mass density of 1.04 x 10 13 g/cm 3 . How- 
ever, this density can be scaled to any desired value n 
by also changing the temperature T so that n/T 3 — 
7.18 x lCT 5 /(0-34360) 3 (MeV-fnLT 3 this insures the value 
of r, see Eqs. 11141 remains the same. Note that our simu- 
lations depend on the electron screening length A only in 
the ratio of A/a. For relativistic electrons, this ratio is in- 
dependent of density. Therefore, the above scaling works 
even with electron screening effects. This is because the 
only length scale in the problem for both electron and 
ion interactions is related to rt -1 / 3 . 

Many run parameters are collected in Table [TTJ We 
evolve the system in time using the simple velocity Ver- 
let algorithm [Tt]] with a time step At — 25 fm/c. We 
use periodic boundary conditions. Our simulation vol- 
ume is large enough so that the box length L = 727.5 fm 
is much larger than the electron screening length A. In- 
deed L/2A = 13.9. The screened potential, for two ions 
separated by a distance L/2, is very small. This helps 
to reduce finite size effects. We include the interactions 
between all particles and do not cutoff the potential at 
large r. We evaluate the interaction between two par- 
ticles as the single interaction with the nearest periodic 
image. We do not include an Ewald sum over further 
periodic images because our box is so large that interac- 
tions with periodic images other than the nearest one are 
very small. 

We start by evolving the system at fixed temperature 
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TABLE II: Simulation Parameters, see text. The energy per 
ion at simulation time t is E(t). 
Parameter 



n 



foof 



E(t = 6x 10 6 fm/c) 
E(t = 151 x 10 6 fm/c) 
(V)/N 
T 



Value 

7.18 x 1(T 
26.17 fm 
328.747886 MeV 
328.747877 MeV 
328.23243 MeV 
0.3436 MeV 



by periodically re-scaling the velocities. We adjust the 
temperature so that approximately half of the system 
remains solid and half liquid. After a simulation time 
of 5 x 10 6 fm/c we switch to evolution at constant en- 
ergy and no longer re-scale the velocities. Thus most of 
our simulation is in the microcanonical ensemble at fixed 
energy and volume. We evolve the system at constant 
energy until the total simulation time is 151 x 10 6 fm/c. 
Energy conservation is excellent. The total energy per 
ion only changed by 3 parts in 10 8 from £ = 6xl 6 fm/c 
to t = 151 x 10 6 fm/c, see Table ILT1 The simulation was 
performed on an accelerated MDGRAPE-2 board [l|| 
and took approximately nine weeks. 



III. RESULTS 

In this section we first test our molecular dynamics 
procedure by simulating a pure system. Then we present 
results for our mixture. A 3456 ion pure system, where 
each ion has the same charge (Z — 29.4) and mass, is 
simulated. One half of the initial configuration is solid 
and the other half is liquid. The system is evolved at 
constant energy for approximately 300000 fm/c. During 
this time, the temperature is expected to evolve to the 
melting temperature because of the release of latent heat, 
as new solid melts or forms. Near the end of the simu- 
lation, we evaluate the temperature as 2/3 of the kinetic 
energy per ion and from this we determine T. We find 
r = 176.1 ± 0.7. The ±0.7 error is statistical only and 
does not include possible errors from finite size or non- 
equilibrium effects. Our result is in good agreement with 
the known T = 175 melting point of the OCP [H]. This 
shows that our molecular dynamics procedure can accu- 
rately describe crystallization, at least for a pure system. 

Next, we calculate the latent heat by determining the 
potential energy difference of 3456 ion pure liquid and 
pure solid configurations. The potential energy difference 
is equal to the latent heat if one assumes the difference in 
density between the phases is small. We find the poten- 
tial energy difference per ion is 0.758 ± 0.002T^/, where 
Tm is the melting temperature. Again, the 0.002 error 
is statistical only and does not include finite size effects. 
Our result is in reasonable agreement with the potential 
energy difference for the OCP of 0.7789T M 13]. Our 
slightly lower melting temperature and latent heat may 
reflect screening length effects in a Yukawa fluid com- 



pared to the OCP [l4J. Alternatively, our slightly lower 
latent heat may reflect finite size effects for a 3456 ion 
system. This latent heat is probably not an important 
heat source compared to the larger energy released from 
nuclear reactions [16j . 

We go on to present results for our mixture with 27648 
ions. The potential energy per ion (V)/N slowly de- 
creases with simulation time until t w 70 x 10 6 fm/c. This 
decrease may be associated with the change in composi- 
tion of the two phases, see below. Next small fluctuations 
are observed in (V)/N for later times that appear to be 
associated with fluctuations in the amount of solid phase 
present in the simulation. The potential energy averaged 
over the last 20 x 10 6 fm/c is given in Table ILT1 The tem- 
perature is evaluated as 2 /3 of the kinetic energy per ion 
and we find T = 0.3436 MeV. 

The parameter T, Eq. [TJ can be evaluated for a mixture 
of ions. For a single ion of charge Zi, the ion sphere radius 
di is the radius of a sphere that contains Zi electrons, 



3Z, 



^Pch 



1/3 



(3) 



with p c h the electron density (or ion charge density). 
Therefore Ti for this ion is, Ti — Zfe 2 /(diT) and av- 
eraging this over a distribution of ions yields T for the 
mixture, 



r = 



(Z 5 / 3 ) 
T 



471-/9, 



'<::/! 



1/3 



(4) 



Note that for a pure system, this equation reduces to Eq. 
□ Table Hn] gives values for (Z 5 / 3 ) and T. The value we 
find for our mixture T = 247 is higher than that for a pure 
OCP (r = 175). This suggests that all of the impurities 
in our crystal phase have somewhat lowered its melting 
temperature. However, see the discussion below about 
chemical separation. 

The configuration of the 27648 ions at the end of the 
simulation is shown in Fig. [3] The solid phase is visible in 
the upper half of the simulation volume where the crystal 
planes are clearly evident. The first interface between 
solid and liquid is just below the center of the box and 
the second interface is near the top of the box. Thus the 
liquid phase extends from the bottom to the top of the 
box because of periodic boundary conditions. Figure [4] 
shows the final configuration of the 832 oxygen ions Z = 
8. The oxygen ions are clearly not distributed uniformly. 
Comparing Fig. [4] to Fig. [3] we see that the oxygen is 
greatly depleted in the solid and enriched in the liquid 
phase. This directly demonstrates phase separation and 
shows that the composition of the liquid is different from 
that of the solid. 

To further explore composition differences, we divide 
the ions into 15 groups according to their z coordinates. 
The first group includes z values from to L/15, etc. The 
average charge of all of the ions in each group is plotted in 
Fig. 03 Groups 1-5 have a relatively small (Z) near (Z) w 
28 while groups 8-13 have a large (Z) r* 30.5. After 
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FIG. 3: (Color on line) Configuration of the 27648 ions at the 
end of the simulation. The crystal planes of the solid phase 
are visible in the upper half of the figure. The lower half of 
the figure shows a liquid phase. The simulation volume is a 
cube 727.5 fm on a side. 



FIG. 5: Average ion charge (Z) in each of 15 sub- volumes. 
Group 1 is at the bottom and group 15 is at the top of the 
simulation volume. 




FIG. 4: (Color on line) Configuration of the 832 oxygen ions 
at the end of the simulation. Oxygen is depleted in the solid 
phase, compare with Fig. [3] 



comparing Fig. [3] with Fig. we somewhat arbitrarily 
identify groups 1-5 as containing liquid phase, groups 8- 
13 solid phase and groups 6-7 and 14-15 as containing 
the two interfaces. See Table Hill 

The composition of the liquid (groups 1-5) and solid 
(groups 8-13) are plotted in Fig. [2l note the log scale, 
and listed in Table [U The compositions of the liquid and 
solid are very different. Chemical elements with Z < 20 
are greatly depleted in the solid phase, while most high Z 
elements are enhanced in the solid phase. Figure [6] plots 
the ratio of the composition in the solid to that in the 
liquid phase for different simulation times. This ratio, 
at t = 151 x 10 6 fm/c, is approximately linear in Z for 



15 < Z < 36. This suggest the affinity of a given element 
for the solid decreases as Z decreases from that of the 
dominant crystal species Z — 34. Elements with even 
smaller Z < 15, while still greatly depleted in the solid, 
do not follow this linear trend. Perhaps very small Z 
ions can occupy interstitial sites in the solid in addition 
to replacing higher Z ions at normal lattice sites. This 
could enhance their concentration in the solid. 

Finally the highest charge ions Z = 47 are, in fact, 
depleted in the solid. This goes against the general rule 
that the solid is enriched in high Z ions. Note that there 
are only a few Z = 47 ions in the simulation. Therefore 
statistical errors could be large. Perhaps this enhance- 
ment of Z = 47 in the liquid is a non-equilibrium effect 
and could go away with further time evolution. However 
we note that for Z = 47 the ratio of solid concentration 
to that in the liquid has been decreasing with simula- 
tion time. Therefore it may be unlikely for the ratio to 
change direction and finally increase with further simu- 
lation time. Instead, the reduction in concentration of 
the solid may be because Z = 47 is a much larger charge 
than the dominant Z — 34 of the crystal lattice. This 
large charge may fit poorly into the existing lattice and 
so the ions may move, instead, into the liquid phase. 

We now address the important question of a further 
time dependence of the composition and if our simulation 
has reached thermodynamic equilibrium. In Fig. [6] we 
plot the ratio of the composition of the solid to that in the 
liquid for different simulation times t. This ratio starts 
at one and decreases, at small Z, with increasing time. 
Comparing the ratio at t = 113 x 10 6 fm/c with that at 
151 x 10 6 fm/c reveals a small but perhaps systematic 
difference. However, we caution that this figure is based 
on our somewhat arbitrary choices of liquid and solid 
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regions at different times. If this difference with time is 
real it may suggest that the composition will continue 
to evolve very slowly for even larger simulation times. 
This is an important open question. In the future we 
will present results for longer simulation times and for 
simulations that start with very different compositions 
for the liquid and solid. Nevertheless, we believe the 
ratio in Fig. [5] clearly shows that the liquid and solid are 
expected to have very different compositions. 



TABLE III: Properties of the original mixture and of the final 
liquid and solid phases after a simulation time of 151 x 10 6 
fm/c, see text. The impurity parameter Q gives the mean 
square dispersion in charge, see Eq. [5] p c h is the ion charge 
density and pb is the baryon density. 

Parameter Mixture Liquid Solid 
{Z) 29.30 28.04 30.48 

Q = {AZf 38.9 52.7 22.3 

(Z 5/3 ) 285.8 269.0 301.5 

(A) 87.62 83.8 91.2 

pch (fm - ' 
Pb (fm~ 3 ) 

T 247 233 261 
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FIG. 6: (Color on line) Ratio of composition in the solid 
phase to that in the liquid phase versus atomic number Z, for 
simulation times t of 10 x 10 6 fm/c (dotted circles) to 151 x 10 6 
fm/c (solid downward pointing red triangles). 



Finally, we discuss the charge and mass densities of 
the two phases, see Table Mil Within small statistical 
errors, we find that the charge density of the liquid is 
equal to that of the solid. This implies that the number 
density of ions is larger in the liquid phase because the 
average ion charge (Z) is larger in the solid phase. This 
equality of charge densities is expected in order to can- 
cel the electron charge density. We find that the average 
mass number (A) is lower in the liquid than in the solid 
phases. Finally the baryon density of the liquid is slightly 
smaller than that of the solid phase. Note that this small 
difference in density may have a significant statistical er- 
ror and may be sensitive to the original distribution of 
Z and A that we use [l6[ • Neutron stars have very large 
gravitational fields. Therefore, chemical separation fol- 
lowed by the sinking of the denser phase, can provide a 
significant source of heating. Although we find only a 
small density difference, this should be checked in future 
work involving different initial compositions. Table IIIII 
also lists values of (Z 5 / 3 ) and T, see Eq. [4] for the differ- 
ent phases. Because (Z 5 ^ 3 ) is smaller in the liquid phase 
we find that T is about 10% smaller in the liquid phase 
compared to that in the solid phase. 



IV. DISCUSSION AND CONCLUSIONS 

How will chemical separation change the structure of 
a neutron star? Consider a steady state situation where 
matter accretes onto a thin ocean while ocean material 
crystallizes to form new neutron star crust. We assume 
the mass of the crust is much larger than that of the 
ocean. In steady state, the rate of crystallization is equal 
to the accretion rate. Furthermore, let us assume the 
composition of the crust is uniform. Steady state equi- 
librium than requires the composition of the crust to be 
equal to that of the accreting material. 

However, the composition of the thin ocean must be- 
come significantly enriched in light elements so that this 
liquid can be in thermodynamic equilibrium with the 
solid crust. Note that the ocean became enriched in light 
elements because the first material to crystallize was de- 
pleted in light elements. Furthermore, this initial change 
in composition of the crystallized material will not no- 
ticeably change the net composition of the crust because 
the crust is assumed to be much more massive than the 
ocean. 

We find a significant enrichment of oxygen Z — 8 in 
the liquid. Our original composition did not include any 
carbon Z = 6 because Gupta et al. [16| found the car- 
bon was burned to oxygen. However if this incorrect 
and carbon is present, it should also be enriched in the 
liquid because it has a similar atomic number to oxy- 
gen. Therefore carbon could be significantly enriched in 
the ocean compared to its concentration in either the 
accreting material or in the crust. Alternatively, carbon 
may burn, either stably or unstably, before it reaches this 
phase transition region. In this case, because there is no 
carbon remaining, it will not be enriched in the liquid. 

Very energetic type I X-ray bursts known as super- 
bursts [l9|, [2(| are thought to involve unstable carbon 
burning. Cumming and Bildsten argue that the mass 
fraction of carbon must be large, X12 ~ 0.05 — 0.10, in 
order for carbon to burn explosively [1, |j| . Chemical sep- 
aration, which we find upon crystallization, could possi- 
bly change carbon concentrations. In addition, chemical 
separation could change the thermal conductivity of the 
crust. This will be discussed in later work, and could 
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impact how the ignition temperature is reached for su- 
perbursts. In addition, the release of latent heat and or 
gravitational potential energy could change the tempera- 
ture profile of the star. However, the small latent heat of 
a pure system, that we found at the beginning of section 
IIII| and the small density difference between our liquid 
and solid phases suggests that both of these heat sources 
may be small. 

We find a lower melting temperature for our mixture 
compared to that for a one component plasma. Table ITTTl 
lists r = 233 for our liquid phase, compared to a pure 
one component plasma, that melts near T = 175. Pre- 
sumably this is due to the large range of charges Z that 
are present in our liquid phase. This change in melting 
point could significantly increase the thickness of the liq- 
uid ocean in accreting neutron stars. If the melting point 
does occur at T = 233, this implies that for accreting neu- 
tron stars with typical crust temperature ~ 5 x 10 8 K that 
the density at which crystallization occurs is, 

p = 2.1 X 10 10 g cm- 3 (T/5 x 10 8 K) 3 (r/233) 3 , (5) 

for the (Z), (Z 5/3 ) and {A) values in table III. A rather 
high density of 2.1 x 10 10 g/cm 3 for crystallization may 
be an order of magnitude higher than the density where 
12 C fuses. Note that the frequency drifts of oscillations 
observed during X-ray bursts may be a way to test the 
depth of the crust/ ocean interface [2lj . 

However, we caution that the melting point could 
change if our simulation is not fully in thermodynamic 
equilibrium. The phase diagram for multi-component 
systems can be very complicated. For example, an addi- 
tional new solid phase could form with a lower Z com- 
position after most of the high Z ions have solidified. 
Therefore it is important to study further the melting 
point of these complex mixtures. In future work we will 
also study phase separation for superburst ashes. 

Our initial composition was not determined in a way 
that is consistent with chemical separation. Gupta et 
al. [l6( calculated how electron capture and light parti- 
cle reactions change the composition of rp process ash as 
it is compressed to higher densities. However, they as- 
sumed the composition does not change upon crystalliza- 
tion. We now find the composition changes significantly. 
Therefore, one should recalculate electron capture and 
light particle reactions consistently with chemical sepa- 
ration. We have calculated results for only one initial 
composition. We expect our general result, that the liq- 
uid is greatly enriched in low Z elements, to hold for 
a variety of different compositions. Nevertheless, it is 
important to study chemical separation for other compo- 
sitions. 

Itoh and Kohyama find the thermal conductivity of 
an impure crystal to be proportional to l/Q where the 
impurity parameter Q is the square of the dispersion in 



the ion charges [22], 

Q = (AZ) 2 = (Z 2 ) - (Z) 2 . (6) 
We find that chemical separation reduces Q from 38.9 in 
the original mixture to 22.3 in the solid phase, see Ta- 
ble [Hll This is because the solid contains far fewer low 
Z ions. Therefore, chemical separation may significantly 
change the thermal conductivity of the crust. Note that 
Q for the liquid phase is also greatly changed. In future 
work we will present molecular dynamics simulation re- 
sults for the static structure factor of both the liquid and 
solid phases and calculations of the thermal conductivity. 

The assumptions, that the composition of the crust is 
uniform and that the system is in steady state equilib- 
rium, are likely to be oversimplified. Instead the crys- 
tallization rate and composition may be time dependent. 
Chemical separation could lead to the formation of lay- 
ers in the crust. There may be bands of high Z mate- 
rial above or below bands of low Z material. This will 
increase the complexity of the crust and it will likely im- 
pact many crust properties. For example, these layers 
could decrease the net thermal conductivity and change 
the temperature profile. Alternatively, if the layers are 
position dependent and dynamically stable, they could 
change the mass quadruple moment of the star and en- 
hance the radiation of continuous gravitational waves. 
The possibility of layers should be studied in future work. 

In conclusion, nucleosynthesis on the surface of accret- 
ing neutron stars likely produces a range of chemical el- 
ements. We have performed molecular dynamics simula- 
tions of crystallization to see how this complex material 
forms new neutron star crust. We find chemical sep- 
aration, with the liquid ocean phase greatly enriched in 
low atomic number elements compared to the solid crust. 
This change in composition can change many crust prop- 
erties such as the thermal conductivity or shear modulus. 
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